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ABSTRACT 

£SJ . The purpose of this communication is to discuss the sim- 
I> ulation of a free surface compressible flow between two fluids, 
typically air and water. We use a two fluid model with the same 
velocity, pressure and temperature for both phases. In such a 
numerical model, the free surface becomes a thin three dimen- 
sional zone. The present method has at least three advantages: 



( i) the free-surface treatment is completely implicit; ( ii) it can 



qq naturally handle wave breaking and other topological changes in 
the flow; ( Hi) one can easily vary the Equation of States (EOS) of 
each fluid (in principle, one can even consider tabulated EOS). 

• Moreover, our model is unconditionally hyperbolic for reason- 
able EOS. 

??: 

Introduction 

One of the challenges in Computational Fluid Dynamics 
(CFD) is to determine efforts exerted by waves on coastal struc- 
tures. Such flows can be quite complicated and in particular 
when the sea is rough, wave breaking can lead to flows that can- 
not be described by basic models like e.g. free surface Euler or 
Navier-Stokes equations. In a free surface model, the boundary 
between the gas (air) and the liquid (water) is a surface. The 
liquid flow is assumed to be incompressible, while the gas is rep- 
resented by a media, above the liquid, where the pressure is con- 
stant (the atmospheric pressure in general). Such a description 
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is known to be valid as far as propagation in the open sea of 
waves with moderate amplitude is concerned. Clearly it is not 
satisfactory when waves either break or hit coastal structures like 
offshore platforms, jetties, piers, breakwaters etc. . . . 

In this work, our goal is to investigate a two-fluid model for 
this kind of problem. It belongs to the family of averaged mod- 
els, that is although the two fluids considered are not miscible, 
there exists a length scale £ in order that each averaging volume 
(of size £ 3 ) contain representative samples of each of the fluids. 
Once the averaging process is done, it is assumed that the two 
fluids share, locally, the same pressure, temperature and veloc- 
ity. Such models are called homogeneous models in the litera- 
ture. They can be seen as limiting case of more general two-fluid 
models where the two fluids could have different temperatures 
and velocities Q. 

The influence of the presence of air in wave impacts is a dif- 
ficult topic. While it is usually thought that the presence of air 
softens the impact pressures, recent results by Bullock et al. 
show that the cushioning effect due to aeration via the increased 
compressibility of the air- water mixture is not necessarily a dom- 
inant effect. First of all, air may become trapped or entrained 
in the water in different ways, for example as a single bubble 
trapped against a wall, or as a column or cloud of small bubbles. 
In addition, it is not clear which quantity is the most appropri- 
ate to measure impacts. For example some researchers pay more 
attention to the pressure impulse than to pressure peaks. The 
pressure impulse is defined as the integral of pressure over the 
short duration of impact. Bagnold 0, for example, noticed that 
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the maximum pressure and impact duration differed from one 
identical wave impact to the next, even in carefully controlled 
laboratory experiments, while the pressure impulse appears to 
be more repeatable. For sure, the simple one-fluid models which 
are commonly used for examining the peak impacts are no longer 
appropriate in the presence of air. There are few studies dealing 
with two-fluid models. An exception is the work by Peregine and 
his collaborators. Wood, Peregrine & Bruce |4] used the pres- 
sure impulse approach to model a trapped air pocket. Peregrine 
& Thais examined the effect of entrained air on a particu- 
lar kind of violent water wave impact by considering a filling 
flow. Bullock et al. [6] found pressure reductions when compar- 
ing wave impact between fresh and salt water where, due to the 
different properties of the bubbles in the two fluids, the aeration 
levels are much higher in salt water than in fresh. H. Bredmose 
recently performed numerical experiments on a two-fluid system 
which is quite similar to the one we will use below. He developed 
a finite volume solver for aerated flows named Flair Q. 



Mathematical model 

In this section we present the equations which govern the 
motion of two phase mixtures in a computational domain £1. 
First of all, we need to introduce the notation which will be 
used throughout this article. We use superscripts ± to denote 
any quantity which is related to liquid and gas respectively. For 
example, a + and a~ denote the volume fraction of liquid and 
gas and obviously satisfy the condition a + + a~ = 1. Then, we 
have the following classical quantities: p ± , u, p, e, E, g which 
denote the density of each phase, the velocity field vector, the 
pressure, the internal & total energy and the acceleration due to 
gravity correspondingly. 

Conservation of mass (one equation for each phase), mo- 
mentum and energy lead to the four following equations: 



ance laws 



(a ± p ± ) / +V-(a ± p ± w) 
(pff)f + V- (pu®u + pl) 
( P E) t + V-(pHu) 



0, 

Pg-u, 



(1) 

(2) 
(3) 



where p := a + p + + oc~ 
specific enthalpy), E = 



p~ (the total density), H := E 



& (the 



\u\ 2 (the total energy). This system 



can be seen as the single energy and infinite drag limit of the 
more conventional six equations model |1|. The above system 
contains five unknowns a ± p ± , u, p and E and only four govern- 
ing equations (Q} - ©. In order to close the system, we need to 
provide the so-called equation of state (EOS) p = p ± (p ± ,e ± ). 
The construction of the EOS will be discussed below. 

It is possible to rewrite these equations as a system of bal- 



dt 



-V-J^(w) =S(w), 



(4) 



where vr(x,t) : M. d x R + R m is the vector of conservative vari- 
ables (in the present study d = 2 or 3 and m = 5), ^(w) is the 
advective flux function and S(w) the source term. 

The conservative variables in the 2D case are defined as fol- 
lows: 

w = (w/)f=i : = (a + p + ,a"p", pu, pv, pE). (5) 

The flux projection on the normal direction n— (ft 1,^2) can be 
expressed in physical and conservative variables 

^ -n — (a + p + %, a~ p~u n ,puu n -\-pn\,pvu n -\-pn2,pHu n ) = 

/ W3fti+W4^2 W3fti+W4^2 \V3ri1-\-W4n2 

k w\— — ,w 2 — — — ,w 3 — — Vpni, 



W\ + W2 W\ + W2 



W4- 



Wl + w 2 



W 3 fti+W 4 ft 2 \ 

+ pn 2 ,(W5+P) ■ (6) 

W\ + W2 / 



where u n := u-n = un\+ vri2 is the velocity projection on the 



<9(J^-H)(w) 
dw 



can 



normal direction n. The jacobian matrix A n (w) := 
be easily computed. This matrix has three distinct eigenvalues: 

Ai = U n — Cy, ^2,3,4 = U n , X$ = U n J r Cy, 

where c s is the sound speed in the mixture. Its expression can 
be found in (8). One can conclude that the system (Q]) - 0]) is 
hyperbolic. This hyperbolicity represents the major advantage of 
this model. The computation of the eigenvectors is trickier but 
can still be performed analytically. We do not give here the final 
expressions since they are cumbersome. 

Equation of state 

In the present work we assume that the light fluid is de- 
scribed by an ideal gas type law 



p =(y-i)p e , 



(7) 



while the heavy fluid is modeled by Tait's law. In the literature 
Tait's law is sometimes called the stiffened gas law l9l[TQli: 



-7T = (^-l)p + ^, 



7To 



^P+' 



(8) 
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where the quantities y, cf, jY , 7Tq are constants. For example, 
pure water is well described when we take JV = 7 and 7Zq — 
2.1 x 10 9 Pa. 

Remark 1 . In practice, the constants cf can be calculated after 
simple algebraic manipulations of equations (0), © and match- 
ing with experimental values at normal conditions: 

_ _ Po + _ ^ Po + Tip 

Cv ~ (y-l)p -r ' Cv ~ (^-l)^p+T ' 

The sound velocities in each phase are given by the follow- 
ing formulas: 

(c;) 2 = r f, ict? = ^±^. (9) 

In order to construct an equation of state for the mixture, 
we make the additional assumption that the two phases are in 
thermodynamic equilibrium: 

p + =p~, T + = T~. (10) 

Below, values of the common pressure and common temperature 
will be denoted by p and T respectively. The technical details 
can be found in Chapter 3, ifTTIl . 

Finite volume scheme on unstructured meshes 

Finite volume methods are a class of discretization schemes 
that have proven highly successful in solving numerically a wide 
class of conservation law systems. These systems often come 
from compressible fluid dynamics. When compared to other dis- 
cretization methods such as finite elements or finite differences, 
the primary interests of finite volume methods are robustness, 
applicability on very general unstructured meshes, and the in- 
trinsic local conservation properties. Hence, with this type of 
discretization, we conserve "exactly" the mass, momentum and 
total energ)Q. 

In order to solve numerically the system of balance laws (Q} 
- © we use ©. The system © should be provided with initial 
condition 

wM)=woM (li) 

and appropriate boundary conditions. 



n K L 




Figure 1 . An example of control volume K with barycenter O. The nor- 
mal pointing from K to L is denoted by Hkl- 

The computational domain £1 C R d is triangulated into a set 
of non overlapping control volumes that completely cover the 
domain. Let 2F denote a tesselation of the domain £1 with control 
volume K such that 

U Ke *rK = Cl, K:=KUdK. 

For two distinct control volumes K and L in the intersection 
is either an edge (2D) or face (3D) with oriented normal Hkl or 
else a set of measure at most d — 2 (in 2D it is just a vertex, in 
3D it can also be a segment, for example). We need to introduce 
the following notation for the neighbourhood of K: 

JV(K) :={Le ST :area(^ HL) ^0}, 

a set of all control volumes L which share a face (or an edge in 
2D) with the given volume K. In this article, we denote by vol(-) 
and area(-) the d and d — 1 dimensional measured respectively. 

The choice of control volume tesselation is flexible in the fi- 
nite volume method. In the present study we retained a so-called 
cell-centered approach, which means that degrees of freedom are 
associated to cell bary centers. 

The first steps in Finite Volume (FV) methods are classi- 
cal. We start by integrating equation © on the control volume K 
(see Figure [T] for illustration) and we apply Gauss-Ostrogradsky 
theorem for advective fluxes. Then, in each control volume, an 
integral conservation law statement is imposed: 

4 / wd£}+ / ^(vf)-n KL dG= [ S(w)dCl. (12) 
at Jk JdK Jk 



lr This statement is true in the absence of source terms and appropriate bound- 
ary conditions. 



2 In other words, in 3D the notation area(-) and vol(-) are very natural and 
mean area and volume respectively, while in 2D they refer to the area and the 
length. 
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Physically an integral conservation law asserts that the rate of 
change of the total amount of a quantity (for example: mass, 
momentum, total energy, etc) with density w in a fixed control 
volume K is balanced by the flux & of the quantity through the 
boundary dK and the production of this quantity y inside the 
control volume. 

The next step consists in introducing the so-called control 
volume cell average for each K G 2? 



that the solution is continuous through an interface) reduces 
to the true flux of the same state, i.e. 

<E>(w,w;H) = • n) (w) . 

After introducing the cell averages wk and numerical fluxes 
into (fT2l) . the integral conservation law statement becomes 



vol(X) 

After the averaging step, the finite volume method can be in- 
terpreted as producing a system of evolution equations for cell 
averages, since 



4 / w(x,t)d& = vol(K) 
at Jk 



d\VK 

dt 



Godunov was the first [ 12 ] who pursued and applied these ideas 
to the discretization of the gas dynamics equations. 

However, the averaging process implies piecewise constant 
solution representation in each control volume with value equal 
to the cell average. The use of such representation makes 
the numerical solution multivalued at control volume interfaces. 
Thereby the calculation of the fluxes Jsk^( w ) ' ^kl da at these 
interfaces is ambiguous. The next fundamental aspect of finite 
volume methods is the idea of substituting the true flux at inter- 
faces by a numerical flux function 

(^( w ) ' n ) IdKndL <— *(w*,w L ;#izL) : M m x R m ^ R m , 

a Lipschitz continuous function of the two interface states 
and Wl. The heart of the matter in finite volume method consists 
in the choice of the numerical flux function <£. In general this 
function is calculated as an exact or even better approximate local 
solution of the Riemann problem posed at these interfaces. In the 
present study we decided to choose the numerical flux function 
according to FVCF scheme extensively described in |[T3l . 
The numerical flux is assumed to satisfy the properties: 

Conservation. This property ensures that fluxes from adjacent 
control volumes sharing an interface exactly cancel when 
summed. This is achieved if the numerical flux function sat- 
isfies the identity 

<Z>(w K ,w L ',n KL ) = -®(w L ,w K ;n LK ). 

Consistency. The consistency is obtained when the numerical 
flux with identical state arguments (in other words it means 



dWK 

dt 



I 



area(Ln/r) 
\ol(K) 



^(wk^l-Hkl) = 



w)L s(w> 



\ol(K) 



d£l. 



We denote by Sk the approximation of the following quantity 
vol |^ J K S(w) d£l. Thus, the following system of ordinary dif- 
ferential equations (ODE) is called a semi-discrete finite volume 
method: 



d\VK 

dt 



vol(K) 



(13) 

The initial condition for this system is given by projecting (ITTb 
onto the space of piecewise constant functions 



w*(0) 



W* woW 



vol(^) 



da. 



This system of ODE (IT3l) should also be discretized. There 
is a variety of explicit and implicit time integration methods. We 
chose the following third order four- stage SSP-RK(3,4) scheme 
CHH with CFL = 2: 





= «W + 






= «« + 






- -11 w 








.(»+!) 







Sign matrix computation 

In the context of the FVCF scheme (see 1 13 ] for more de- 
tails), we need to compute the so-called sign matrix which is 
defined in the following way 

U n := sign(A„) =7?sign(A)L, 



4 



Copyright © 2008 by ASME 



where R, L are matrices composed of right and left eigenvectors 
correspondingly, and A = diag(Ai , . . . , A m ) is the diagonal matrix 
of eigenvalues of the Jacobian. 

This definition gives the first "direct" method of sign matrix 
computation. Since the advection operator is relatively simple, 
after a few tricks, we can succeed in computing analytically the 
matrices R and L. For more complicated two-phase models it is 
almost impossible to perform this computation in closed analyti- 
cal form. In this case, one has to apply numerical techniques for 
eigensystem computations. It turns out to be costly and not very 
accurate. In the present work we use physical information about 
the model in the numerical computations. 

There is another way which is less expensive. The main idea 
is to construct a kind of interpolation polynomial which takes the 
following values 

P(u n ±c s ) = sign(w„±c iS ), P{u n ) = sign(w„). 

These three conditions allow us to construct a second degree in- 
terpolation polynomial. Obviously, when P{X) is evaluated at 
A = A n we obtain the sign matrix U n as a result. The construc- 
tion of the Lagrange interpolation polynomial P(A) is simple. 

In our research code we have implemented both methods. 
Our experience shows that the interpolation method is quicker 
and gives correct results in most test cases. However, when we 
approach pure phase states, it shows a rather bad numerical be- 
haviour. It can lead to instabilities and diminish overall code 
robustness. Thus, whenever possible we suggest to use the com- 
putation of the Jacobian eigenstructure. 

Second order scheme 

If we analyze the above scheme, we understand that in fact, 
we have only one degree of freedom per data storage location. 
Hence, it seems that we can expect to be first order accurate 
at most. In the numerical community first order schemes are 
generally considered to be too inaccurate for most quantitative 
calculations. Of course, we can always make the mesh spacing 
extremely small but it cannot be a solution since it makes the 
scheme inefficient. From the theoretical point of view the situa- 
tion is even worse since an ff{h2) Li-norm error bound for the 
monotone and E-flux schemes [|T6l is known to be sharp [17], 
although an G{H) solution error is routinely observed in numer- 
ical experiments. On the other hand, Godunov has shown [12] 
that all linear schemes that preserve solution monotonicity are at 
most first order accurate. This rather negative result suggests that 
a higher order accurate scheme has to be essentially nonlinear in 
order to attain simultaneously a monotone resolution of discon- 
tinuities and high order accuracy in continuous regions. 

A significant breakthrough in the generalization of finite vol- 
ume methods to higher order accuracy is due to Kolgan [18, 19 ] 



and van Leer [20]. They proposed a kind of post- treatment pro- 
cedure currently known as solution reconstruction or MUSClE 
scheme. In the above papers the authors used linear reconstruc- 
tion (it will be retained in this study as well) but this method was 
already extended to quadratic approximations in each cell. 

In this paper we briefly describe the construction and prac- 
tical implementation of a second-order nonlinear scheme on un- 
structured (possibly highly distorted) meshes. The main idea is 
to find our solution as a piecewise affine function on each cell. 
This kind of linear reconstruction operators on simplicial con- 
trol volumes often exploit the fact that the cell average is also 
a pointwise value of any valid (conservative) linear reconstruc- 
tion evaluated at the gravity center of a simplex. This reduces 
the linear reconstruction problem to that of gradient estimation 
at cell centers given cell averaged data. In this case, we express 
the reconstruction in the form 

w*(5) = + (Vw) z -(x-xo), K e ST , (14) 

where wk is the cell averaged value given by the finite volume 
method, (Vw)^ is the solution gradient estimate (to be deter- 
mined) on the cell K,xeK and the point xo is chosen to be the 
gravity center for the simplex K. 

It is very important to note that with this type of representa- 
tion (fT4l) we remain absolutely conservative, i.e. 

! F v / yvK(x)da = w K 

vo\(K) Jk 

due to the choice of the point xq. This point is crucial for fi- 
nite volumes because of intrinsic conservative properties of this 
method. 

In our numerical code we implemented two common tech- 
niques of gradient reconstruction: Green-Gauss integration and 
least squares methods. In this paper we describe only the least 
squares reconstruction method. The Barth-Jespersen limiter [21] 
is incorporated to obtain non-oscillatory resolution of disconti- 
nuities and steep gradients. We refer to [fTTTl for more details. 

Least-squares gradient reconstruction method 

In this section we consider a triangle control volume K with 
three adjacent neighbors T\, and T^. Their bary centers are 
denoted by O(xo), 0\ (x\), 02(^2) and 0^(x^) respectively. In the 
following we denote by w; the solution value at gravity centers 
Of. 

v?i \=w(xi), w := w(x ). 



3 MUSCL stands for Monotone Upstream-Centered Scheme for Conservation 
Laws. 
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the left both sides of $M by [LiL 2 ] f yields 



2 



. O 



3 



K 



.0, 



G(Vw)* = fc, G=(Z, 



ijh<i,j<l ' 



(Li -Li) (Li-L 2 )\ 



(Lz-Li) (L 2 -L 2 )J 



(19) 



where G is the Gram matrix of vectors 



{L1X2} 



and b 



( "0 ) . The so-called normal equation dT9b is easily solved 
by Cramer's rule to give the following result 



Figure 2. Illustration for least-squares gradient reconstruction. We de- 
pict a triangle control volume with three adjacent neighbors. 



(Vw)i 



J (hiiLi'h-hiiLi'f) 



hih2-l 2 n \ln{L2'f)-hi{Li'f) 



Our purpose here is to estimate Vw = (d x w, d y v?) on the cell 
K. Using Taylor formula, we can write down the three following 
relations: 



wi - 

W2 " 

w 3 - 



w 
•w 
•w 



(Vw)^-(xi-x ) + ^ 2 ), 
(Vw)^-(x 2 -xo) + ^(/z 2 ), 



(Vw)#-(x 3 -x ) 



(15) 
(16) 
(17) 



If we drop higher order terms &(h 2 ), these relations can be 
viewed as a linear system of three equations for two unknown^ 
(d x w,dyw). This situation is due to the fact that the number of 
edges incident to a simplex mesh in R d is greater or equal to d 
thereby producing linear constraint equations (1131) - (1771) which 
will be solved analytically here in a least squares sense. 

First of all, each constraint (IT31) - (IT71) is multiplied by a 
weight C0i G (0, 1) which will be chosen below to account for dis- 
torted meshes. In matrix form our non-square system becomes 



0)i Ax\ CGiAyi 
\ CO2AX2 C02AV2 
co 3 Ax 3 (O x Ay 3j 



(VW)* : 



'&>i(wi — Wq)^ 

0>2(W2 —Wq) 
,C03(W3-Wo)> 



where Axi = xi — xo, Ayt = yt — yo. For further developments it is 
convenient to rewrite our constraints in abstract form 



[Li,L 2 ]-(Vw)* = /. 



(18) 



We use a normal equation technique in order to solve symboli- 
cally this abstract form in a least squares sense. Multiplying on 



4 This simple estimation is done for scalar case only w = (w) . For more gen- 
eral vector problems the numbers of equations and unknowns have to be changed 
depending on the dimension of vector w. 



The form of this solution suggests that the least squares linear 
reconstruction can be efficiently computed without the need for 
storing a non- square matrix. 

Now we have to discuss the choice of weight coefficients 
{o)i}^ =l . The basic idea is to attribute bigger weights to cells 
barycenters closer to the node N under consideration. One of the 
possible choices consists in taking a harmonic mean of respective 
distances r ; = — xn\\. This purely metric argument takes the 
following mathematical form: 



\\Xi-X N \\ 

where k in practice is taken to be one or two (in our numerical 
code we choose k = 1). 



Numerical results: falling water column 

A classical test in violent flows is the dam break problem. 
This problem can be simplified as follows: a water column is 
released at time t = and falls under gravity. In addition, there 
is a step in the bottom. During its fall, the liquid hits this step 
and recirculation is generated behind the step. Then the liquid 
hits the vertical wall and climbs along this wall. The geometry 
and initial condition for this test case are depicted on Figure 
Initially the velocity field is taken to be zero. The unstructured 
triangular grid used in this computation contained about 108000 
control volumes (which in this case are triangles). The results 
of this simulation are presented on Figures [4]— El We emphasize 
here that there is no interface between the liquid and the gas. The 
dark mixture contains mostly liquid (90%) and the light mixture 
contains mostly gas (90%). An interesting quantity is the impact 
pressure along the wall. It is shown in Figure [51 where the max- 
imal pressure on the right wall is plotted as a function of time 
f i-> max( xy) e ! x [o, 1] p 0, y, t) . 



6 



Copyright © 2008 by ASME 



Then, we performed other computations with volume frac- 
tions closer to pure phases. For example, we show some results 
with the gas mixture modelled with a + = 0.05, a~ =0.95. The 
pressure is recorded as well and this result is plotted on Fig- 
ure [TUJ One can see that the peak value is higher and the impact 
is more localised in time. 



1 









a+ =0.1 




a" =0.9 


a+ =0.9 




a- =0.1 






n i 005 



0.3 0.65 0.7 1 

Figure 3. Falling water column test case. Geometry and initial condition. 




(a) t = 0.005 (b) t = 0.06 

Figure 4. Initial condition and the beginning of column dropping process. 




(a) t = 0.1 (b) f = 0.125 

Figure 5. Splash creation due to the interaction with the step. 



(a) t = 0.2 (b) t = 0.225 

Figure 6. The liquid hits the wall. 




(a) t = 0.3 (b) t = 0.4 



Figure 7. The splash is climbing along the wall. 




(a) t = 0.5 (b) t = 0.675 



Figure 8. Turbulent mixing process. 



Conclusions 

In this article we presented a simple mathematical model 
for simulating water wave impacts. The preliminary results are 
encouraging and the validation of this approach will be the sub- 
ject of future work. Namely, we are going to perform qualitative 
and quantitative comparisons with the more general six equations 
model Q. 

We also presented an efficient numerical approach for dis- 
cretizing the governing equations. It is a second order finite 
volume scheme on unstructured meshes. This method was im- 
plemented in our research code. By construction, our code has 
excellent mass, momentum and energy conservation properties. 
Numerical tests presented in [1 1 j partially validate the method. 

We also plan to carry out a parametric study with our solver. 
The influence of aeration, gas properties and other factors on the 
impact pressures is very important for industrial applications. 
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Maximal pressure on the right wall 



' 0.1 0.2 0.3 0.4 0.5 0.6 0.7 

t, time 



Figure 9. Maximal pressure on the right wall. Heavy gas case. 

Maximal pressure on the right wall 



' 0.1 0.2 0.3 0.4 0.5 0.6 0.7 

t, time 



Figure 1 0. Maximal pressure on the right wall as a function of time. Light 
gas. 
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